GREEN'S FUNCTION FOR SYMMETRIC LOADING OF AN ELASTIC 
SPHERE WITH APPLICATION TO CONTACT PROBLEMS 



A. S. TITOVICH AND A. N. NORRIS 

Abstract. A compact form for the static Green's function for symmetric loading of an elastic 
sphere is derived. The expression captures the singularity in closed form using standard func- 
tions and quickly convergent series. Applications to problems involving contact between elastic 
spheres are discussed. An exact solution for a point load on a sphere is presented and subse- 
quently generalized for distributed loads. Examples for constant and Hertzian-type distributed 
loads are provided, where the latter is also compared to the Hertz contact theory for identical 
spheres. The results show that the form of the loading assumed in Hertz contact theory is valid 
for contact angles up to about 10 degrees. For larger angles, the actual displacement is smaller 
and the contact surface is no longer flat. 



1. Introduction 

Contact between spheres has intrigued researchers for more than a century, and still no simple 
closed-form analytical solution exists. One of the first and most important developments in the 
field, due to Heinrich Hertz [1] in 1881, is an approximate solution for the normal, frictionless 
contact of linear elastic spheres. The major assumption in Hertz's model is that the contact 
area was small compared to the radii of curvature, which has served as a useful engineering 
approximation in many applications. Ever since then many have tried to relax this assumption 
while maintaining a compact, workable solution. The Green's function for symmetric loading on 
a sphere provides the means to find the exact response for arbitrary loading, a first step towards 
improving on Hertz' classic solution. Existing forms of the Green's function are however not 
suitable for fast and ready computation, either due to slow convergence of series or analytically 
cumbersome expressions. The goal of the present paper is to provide an alternative form of the 
Green's function suitable for fast computation of solutions under arbitrary loading. 

Sternberg and Rosenthal [2] present an in-depth study of the nature of the singularities on 
elastic sphere loaded by two opposing concentrated point forces. As expected, the dominant 
inverse square singularity in the stress components can be removed by subtraction of an appro- 
priate multiple of Boussinesq's solution for a point load at the surface of a half space. Sternberg 
and Rosenthal showed that the quickly convergent residual field retains a weaker singularity 
of logarithmic form, a result that is also evident in the solution developed here. The singular 
solutions obtained by Sternberg and Rosenthal were extended to arbitrarily oriented point forces 
by Guerrero et al. [3]. Our interest here is in developing an analogous separation of the Green's 
function (circular ring loading). In this regard, a relatively compact form of the Green's function 
for the sphere was derived by Bondareva |1] who used it to solve the problem of the weighted 
sphere. In [5], Bondareva formulates an example with a sphere contacting a rigid surface. This 
has been used to solve the rebound of a sphere from a surface [6]. Bondareva's solution starts 
with the known series expansion [7] for the solution of the elasticity problem of a sphere, and 
replaces it with finite integrals of known functions. 
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In this paper we introduce an alternative form for the Green's function for a sphere, comprised 
of analytical functions and a quickly convergent series. No direct integration is required. The 
methodology for determining the analytical functions is motivated by the simple example of a 
point load on a sphere, for which we derive a solution similar in spirit to that of [2], but using 
a fundamentally different approach: partial summation of infinite series as compared with a 
functional ansatz. The present methods allows us to readily generalize the point load solution to 
arbitrary symmetric normal loading. A typical contact problem involves solving a complicated 
integral equation for the contact stress once a displacement is specified. Instead, we will use 
the derived Green's function in the direct sense, solving for the displacements for a given load. 
This is used to check the validity of Hertz contact theory through the assumed form of the stress 
distribution. 

The outline of the paper is as follows. The known series solution for symmetric loading on 
a sphere is reviewed in ^ The proposed method for simplification is first illustrated in fj3] by 
deriving a quickly convergent form of the solution for a point force. The Green's function for 
symmetric loading is then developed in ^ and is illustrated by application to different loadings. 
Conclusions are given in ^ 



2. Series solution 

Consider a solid sphere of radius R, with surface r = R, < 9 < ir, in spherical polar 
coordinates (r, 0, (p). The sphere is linear elastic with shear modulus G and Poisson's ratio u. 
The surface is subject to tractions 

= 0, 0-^0 = 0, a„ = cr{9), for r = i?, < < vr, < < 27r. (1) 

Using the known properties of Legendre functions, see. eqs. (j45p . allows us to express the normal 
stress as 

rr(f)) = ^ 



(y{e) = \Y.^2n + l)a„P„(cos0), (2) 

n=0 

where the Legendre series coefficients are 



j a{(j))Pn{cos(j)) sin (j)d(j). (3) 

The displacements and tractions f( 

5] 



<7n, 



The displacements and tractions for the sphere can also be represented in series form [8l eq. 



2Gur = [{n - 2 + 4i^)Anr + Bnr~'^] r"P„ (cos 6*), 

n=0 

2Gue = V [n(n + 5 - 4z.)Ar + (n + l)Bnr-'] r" 4^^, 
^-^ ^ ^ nin + 1) 

(4) 

= [[n{n - 1) - 2(1 + u)]An + (n - l)Bnr~^] r"P„ (cos 



n=0 
oo 



^re = V [n[(n - l)(n + 3) + 2(1 + u)]An + {n' - l)i?„r"2] r" 4^^^, 
^-^ ^ ^ nin + 1) 

71=1 ^ ' 
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with Bq = 0, and Bi corresponds to a rigid body translation via 2Gu(0, •) = Bie^- It follows 
from that 

A C7„ (n+l)(2n + l) 

4i?" [n(n - 1) + (2n + 1)(1 + z.)] ' " ' ^ ' 

Bn = [{n - l)(n + 3) + 2(1 + z^)] R^An, n > 2. 

— 1 

Thus, noting that P^(cos6l) = ^PnicosO), we have 

R (2{l-2v) ^ (2n + l\(2{l-u)n'^ + vn-l + 2v 



(6a) 

( 2n^\ \( {-l + 2u)n + 2-v \ „ ^ 
= iG g/" J + (1 + 2^)n + 1 + . J 

Bondareva [4j, using a different representation, replaced the infinite summation of Legendre 
functions by a combination of closed form expressions and an integral, each dependent on v. 
The integral term contains a logarithmic singularity which, together with the complex-valued 
nature of its coefficients, makes its evaluation indirect. Here we propose an alternative form for 
the Green's function in a combination of closed-form expressions and a standard summation of 
Legendre functions that is, by design, quickly convergent. 



3. Point force 

3.1. Exact solution. In order to illustrate the method, we first consider the simpler problem 
of the point force of magnitude F applied at = defined by 

-F ^. 5{e-i^) -F 
2'kR^ smV' 2-kR^ 

where we have used the property -Pn(l) = 1- The difficulty with the infinite summations ([6]) is 
two fold: first, it is not a suitable form to reproduce the singular nature of the Green's function; 
secondly, it does not converge quickly as a function of the truncated value for n. The idea here 
is to replace the summation by closed form expressions plus a summation that is both regular 
and quickly convergent. The fundamental idea behind the present method is to write n^, uq of 
eqs. ([6]) in the form 

M 



— F / \ 
Ur{R. 0) = —— 4(1 - ,.)S{e) + ^ a,{,^)S,{e) + f{e) , (Sa) 

V j=o ^ 
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where the functions S{9) and Sj{6) (j = 1, . . . M), are closed- form expressions, in this case, 

oo 19 
Si9) = Y,Pn{cos9) = -csc-, (9a) 



2 2' 

n=0 



S,{9) = f:^^±l^, j=0,l,..., (9b) 



n=0 

and /(^), g{9) are regular functions of 9 defined by quickly convergent series in n, 

oo oo 

f{9) = CnPnicos 9), g{9) = J2 DnPnicos 9). (10) 

n=0 n=0 

The coefficients oq, ai ... ajv/ are defined so that C„ =0(n-(*^+2)) as n — 7- CX3. This criterion 
uniquely provides the constants oq, ai ... om as solutions of a system of linear equations. 
Similarly, bo, bi . . . bu are uniquely defined by Dn =0{n~^^^^'^'>) as n — )• oo. 

Here we consider the specific case of M = 2. Other values of M could be treated in the 
same manner; however, we will show that M = 2 is adequate for the purpose of improving 
convergence. In this case eq. ([8]) becomes 

Ur{R, 9) [4(1 - y)S{9) + aoSo{9) + aiSi{9) + a2S2{9) + f{9) 

~^ lS~^ f A/-, \ "0 CLi 02 



[ V (4(1 - z.) + ^ + ^ + + Cn)Pn{e) + CoPoie) + c^Pi{e) 

L ^-^ \ n + 1 n n — 1 J 



SttGR L ^ V ' ' n + 1 

n=2 

+ Ml-iy)[Po{0) + Pi{9)) +ao(^Po{9) + ^Pi{9)] +a,Pi{9) 



flla) 



-F d 

-F d bo 



ueiR, 9) - \hoSo{9) + biSi{9) + b2S2{9) + g{9) 



SttGR d9 _ ^ 

n=2 



L Vn + l n n — 1 J 



+ 6o(^'o(^) + ^Pi(^)) +6iPi(e) + I?o^'o(e) + ^iA(^)], (lib) 
where the associated three functions Sj{9), j = 0, 1, 2 are (see the Appendix) 



5o(^) = log(^l + csc-J, (12a) 

Si{9) = -So{9)-2logsm^, (12b) 

9 / 9\ 

S2{9) = Si{9) cos - 2 sin - 1 - sin - . (12c) 



Equations ([8]), (I9ap and (I12p indicate the expected Boussinesq-like singularity as well as 
the weaker log^ singularity first described by Sternberg and Rosenthal [2]. The logarithmic 
singularities in Sj{9), j = 0, 1, 2 can be compared to the potential functions [Di], [D2], and [D3] 
in equation (17) of [2], which provide a logarithmic singularity. In the present notation these 
are, respectively (using capital $s so as not to be confused with the angle (p, and making the 
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substitution 9 ^ n — 6) 

6 9 
^>i(0) = 21ogsin-, $2(6') = + 2 cos 61 log sin-), (13a) 

Q 

$3(6!) = R^{2{1 - 2>cos^ 9)\ogs\n- + cos^ 9 - 3 cos 6* - l). (13b) 

These clearly display the same form of the singularity as in equations ()12p . but are otherwise 
different. 




90 

9 (degrees) 



180 



Figure 1. Convergence of the proposed expression for Ur = ^^^ Ur in eq. (jllj] 
The magnified view at 9 = ^ shows the difference in the value of Ur, which is 
iUr)N=io - {Ur)N=4 = 4.5259 • lO'^ and (C/.)7v=ioo - {Ur)N=io = -1.7059 • IQ-^. 
The small plot is magnified by a factor of 11250. 



Define the first two coefficients of f{9) and g{9) from (jlOp as 
2(1 - 2u 



4(l-z^)-ao, Ci = -4{1 - u) - -ao - ai, 







50, 



(14a) 
(14b) 



The coefficients a„ and 6„ are then found by comparing expression (jllh to the series solution 
in dl]), expanding both expressions for large n, and equating the coefficients of the same order 
terms. Thus, the original assumed form of the solution ([8]) implies 



E 

n=2 



Pni9) 



2n+l\ { {2{l - v)n'^ + vn - I + 2v) 



re2 + (1 + 2u)n +l + v\ n-l J y ((-1 + 2v)n + 2-v) 

\\n+l n n-1 ^"■J^ 



(15) 



n=2 



6 



A. S. TITOVICH AND A. N. NORRIS 



where 



ao = ^{l + u){l- 2u){-l6u^ + 8U + 5), 

ai = -32u^ + 16i/^ + 30u^ - 16u - 1, 

13 1 

02 = 16z^^ - 16z^3 - 5z^2 + J, + 

6o = -(l + z^) (161/2 -121.-1), 
6^ = _16i/3 + 4j,2_^i3j,_4^ 



5 5 

'2 



62 = 81/^ - 6z/2 - -z/ + 



2 2 

The remaining coefficients C„ and D„ are then determined directly from (|15p 
Cn = ^ '^'^ ((6 - ai - 6a2)n + ai 



Ln 



(^(64z/^ - 80z/^ + 22,v - 4)n + {2,2v^ - IQv^ - SOi/^ + IQv + 1)) , (17a) 



I?n = ^i^((6-fei-662)n 



^i^tZl ((321/3 _ 32j,2 _ + 5)n + (iGi/^ _ 41,2 _ ^^^^ _^ 4)^ ^ (i^^) 



where 



Ln = n{r? -l){r? + {l + 2u)n + l + v). (18) 

In summary, the new form of the point force solution is given by the displacements in eq. (jlip 
where the functions and coefficients are given in eqs. (|12p -(|14p. (jl6p -(jl8 p . 

It is useful to note that the present approach, like that of Sternberg and Rosenthal [2], leads 
to a separation of the solution into a part that displays the singular behavior plus an additional 
part that is quickly convergent. Sternberg and Rosenthal's solution was based on an ansatz 
[21 eq. (38) ] motivated by separation of the stress into a Boussinesq term with leading order 
singularity plus a residual field. Our approach is somewhat different, in that we regularize the 
infinite series solution directly by partial summation. The net effect is essentially the same 
as in (21 eq. (38) ] for the point load. However, the present method can be easily adapted to 
the more general case of the Green's function for a circular (ring) load, see 21 The analogous 
generalization of Sternberg and Rosenthal's method is not as obvious. 

3.2. Numerical examples. In the following examples we introduce the integer as the trun- 
cation value of the series in eqs. (jlOp . The Poisson's ratio was taken to be 0.4. Displacements 
have been normalized by the constant coefficient of the series as Ui = -SnGRF-^Ui, where 
i = r,6. Figures [D and [2] show the rate of convergence of the displacements given by eq. (Ilip . 
whereas Figures [3] and H] compare the displacements in eqs. ([6]) with eqs. (Ilip . 

By design, the proposed expression (eq. (Ilip ) converges much faster than the existing expres- 
sion (eq. ([6])) as seen in Figures [3] and |H Looking at the convergence of the proposed expressions 
with the truncation value A^, Figures [1] and [21 we can suggest that the analytic portion of the 
expression alone gives close results. However, it should be noted that one cannot get rid of the 
first two terms in the series for f{9) and g{6) because of their large magnitudes. As far as the 
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general behaviour of the normahzed displacements with 0, we see that they increase asymptoti- 
cally approaching = 0, change sign between 36.7° and 108.7° for Ur (7.27° and 80.83° for Ue), 
and have a minimum at 65.5° for Ur (24.6° for Uq). This is difficult to see in the Figures, but 
due to symmetry of the loading, the displacement Uq must have a value of zero at = 0. 



1.5 ■ 



0.5 





N=4 




N=10 




N=100 







Detail at 9=90 



30 60 90 120 

9 (degrees) 



150 



180 



Figure 2. Convergence of Ug = ^^^ ug given in eq. (Ilip . The magnified view 
at 6 = ^ shows the difference in the value of Ug, which is {Ug)N=io — {Ug)N=4: = 
2.4097 • 10~3 and {Ug)N=ioo - {Ug)N=io = -1.5040 • 10^1 The small plot is 
magnified by a factor of 1200. 





30 60 90 120 150 180 




30 60 90 120 150 18 
6 (degrees) 



Figure 3. Comparison of the convergence of Ur = ~^^^ Ur with the truncation 
value for the existing expression (eq. ([6]), shown in blue) and the expression 
proposed herein (eq. (jlip . shown in red). Left to right: A^=4, 10, 100. 
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9 (degrees) 9 (degrees) 6 (degrees) 



Figure 4. Comparison of the convergence of Ug = p uq with the truncation 
value for the existing expression (eq. ([6]), shown in blue) and the expression 
proposed herein (eq. (Ilip . shown in red). Left to right: N=A, 10, 100. 



4. Green's function 

4.1. A fast convergent form for the Green's function. The surface displacements for 
arbitrary loading may be written, by analogy with the ansatz ([5]) for the point force, and 
generalizing the latter, 

M 



MR,0) = [Ylb,{y)S,{e,4>)+9{0A)^ a{ct>) sin cpd^, (19b) 

where 5(6*, (p) and Sj{e, cp) (j = 1, . . . M), are 

00 

S{0, cl)) = J2 ^"(cos ^)i^n,(cos (j)), (20a) 

n=0 
00 ^ 

Sj{e,4>) = -— -Pn+j {cose) Pn+j (cos (j)), J = 0, 1, ... , (20b) 
n=0 

and f{9, (p), g{6, (p) are regular functions of 9 defined by quickly convergent series in n, 

00 00 

f{e,(t>) = Y,CnPn{cOSe)Pn{cOs4>), g{e ,((>) = ^nPn (cOS 9) PnicOS (j)) . (21) 
n=0 n=0 

The coefficients ao, ai • • • aM are the same as before. The main complication is to find the 
functions (f20l) . Thus, S{6,(f)) follows from (f52l) as 

ficscf sec|K('(cot|tan|)^y 6 >(/), 
"^^^'"^^ ~ I isecf csc|K((tan|cot|)^], 6 < (/>, ^^^^ 
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where K{m) is the complete elhptic integral of the first kind [9j 17.3.1], while eq. (i50|) implies 



5, 



on 



5o( 



+ Si 



-2 log sin I cos |, 6 > (f), 
-2 log cos I sin ^, < (f). 



(23) 



The functions Sj{0, (p) can be determined, but their form is overly complicated, and defeats our 
objective of simplifying the Green's function. We therefore restrict the solution to the use of the 
above two series: S{9,(j)) and Soi{9,(l)). 

We therefore consider the following form of the ansatz (jl9|) using the series S and Sqi of eqs. 
(f22|l and (f23]) . respectively. Substituting them into eqs. ([19]) yields the identities 



Ur{R,e) 



ue{R,\ 



E_ 

''AG 



'0 



4(1 - u)S{e, <P) + aoi{u)Soi{e, <P) + f{9, 4>) ) (7(0) sin <pd(t>, 



R_ r 



n=2 



f2 (4(l-i^) + aoi 



2n + l 
n{n + 1) 



+ Cn )-P„(cos6')P„(cosi 



+ 4(1 - i/)( Po(cos6l)Po(cos (/>) + Pi(cos6l)Pi(cos. 



+ aoi ^-Po(cos6')Po(cos(/)) + -Pi (cos 6^) Pi (cos (fy 
+ CoPo(cos6i)Po(cos(/>) + Ci Pi (cos 61) Pi (cos. 



cr up sm 



R_d_ r 



n=2 



2n -)- 1 

'oil — — — ] +Dn )P„(cos6l)P„(cos. 



n(n + 1) 



+ 601(^^0 (cos 6l)Po (cos (/>) + -Pi(cos6')Pi(cos. 
+ P»o^'o(cos6l)Po(cos</>) + P»iPi(cos6l)Pi(cos<; 



cr((/)) sin 



(24a) 



(24b) 



Once again we define the first two coefficients of f{9, (p) and g{9, 0) as 



= ^ y - 4(1 -u)- aoi, Ci = -4(1 -ly)- ^aoi. 



-Do = -boi, 



Di = "2^01, 



(25a) 
(25b) 
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which aUows us to solve the following expressions for the coefficients aoi and 601 • 

r — 2 I n I o \ — TT"^ P„(cos6l)P„(cos(/)) 

\ n — 1 J \ + + Iv)n + 1 + 1/ J 

(2n + l) 



V (4(1 - z.) + aoi + C„) P„ (cos (cose/.), (26a) 

V nyn + 1) / 



n=2 



,^/2n + l\/ (-l + 2i/)n + 2-z^ \ , , 
V ( — I ( -i — ^ '—^ P„(cos e)Pn(cos ( 



n - 1 / Vn2 + (1 + 2z/)n + 1 + 

n=2 

= £ f^T^VTT^oi + i?„) Pn (cos 0)P„ (cose/.). (26b) 

This is done by expanding equations (I26p for large n and equating same order terms yielding 

aoi = (2z/-l)2, 6oi = 2i/-l. (27) 
Using ([27]), C„ and D n are found directly from equation (I26p (see also (llSp ) 

C„ = - -3- [(1/ - l)(4i/ - l)(4z/ + l)n3 + (8z^^ - llu - l)n^ 

+ (-121/3 _^ + _ 5)j^ - (1/ + l)(2i/ - 1)2] , (28a) 

Dn = - [(z. - l)(4i/ + l)n2 + 2(-z/2 + J, _ i)n - (i/ + l)(2i/ - 1)] . (28b) 



In summary, 



UiiR,0)=— Hi{e,^)a{cl)) sin cl)d(P, i = r,e, (29a) 

00 

/?,(&,(/)) =4(1 - u)S{9,^) + (1 - 2z/)25oi(^,(/.) + J^C„P„(cos0)P„(cos(/>), (29b) 

n=0 

He{e, (f) =— i {2v - l)Sm{e, 0) + DM^os e)P„(cos </.) J (29c) 

^ n=0 ^ 

where the coefficients C,i, are given in ()28p . Note that C^, -Dn =0(n~2) as n — )• 00, ensuring 
rapidly convergent series. The Green's functions of ()29p are generally valid for 9 E [0, vr]. The 
integrands Hi{6,(j)) are smooth and bounded functions of (/> for </> / 0, which is always the case 
if the displacements are evaluated at points outside the region of the loading a{^). However, for 
points under the load, the integration of Hr{0, <j)) involves a logarithmic singularity at (/) = 0. A 
simple means of dealing with this is described next. 



4.1.1. Removing the singularity under the load. The function S{0,(j)) exhibits a logarithmic sin- 
gularity by virtue of the asymptotic behavior 



iC(m) = log : + 0(l-m), m t 1- 

yl — m 



(30) 
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Figure 5. Comparison of the proposed solution (eq. ()35p shown in red) and 
existing series solutions (eq. ([B]) shown in blue) for Ur = ~^^^ Ur under a 
Hertzian-type load distributed up to (p = 10°. Left to right: N=4, 10, 100. 
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Figure 6. Comparison of the proposed solution (eq. ()35p shown in red) and 
existing series solutions (eq. ^ shown in blue) for Ug = ~^^^ ug under a 
Hertzian-type load distributed up to = 10°. Left to right: N=4, 10, 100. 



The integral in (j24ap is evaluated by rewriting eq. ()24ap in the equivalent form 



UriR,e) 



'AG 



4>Q 



{aQi{v)Sm{e,(t>) + f{e,(t>))a{ct>) 



+ A{i-u){s{e,<t>)a{cp)-s{eA)'y{0)) 

+ A{l-v)a{6) Z'^" 5(0,0) sin < < (/>o 



(31) 



where the angle (f)Q defines the domain of the loading, which is normally for contact problems, 
much less that vr. The function S{9,(f)) has the same singularity as S{9,<j)) and has a relatively 
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simple integral. We choose 



2^ CSC I sec I log(cos^ | — cos^ |), 9 > 

1 9 5 1 / ■ 2 <^ -2 

2^: sec 2 CSC S log(sm S — sm 



2^' 



< 



(32) 



The integrand of the first integral in (j3ip is now a smoothly varying function with no singularity, 
and the second integral is, explicitly. 



S{0,(p) sin (j)d4> 



vr sm ■ 



cos I 



log(x^ — cos^ |)(ix 



vr cos 2 



log(x^ — sin^ |)dx 



G(cosf,l) G(sin|,sinf) 



+ 



o<e< 



where 



sm : 



cos ■ 



(33) 



G{x, y)= {{y- x) \og{y -x) + {y + x) log(y + x) - 2{y - x + x log 2x)) . (34) 



vr 



In summary, the solution for Ur with the singularity removed has the following form (see also 
eq. (pa for Hr{9,<i)) and eq. ^ for G{x,y)) 



Hr{e,<p)a{cp)-Hr{e,cp)a{e) 



sm 



+ 



/i(0)}, 



Hr 



--4{i-u)s{e, 



(35) 



/i(0) =4(1 - z/) 



G(cos|,l) G(sinf,sinf) 



+ 



sm : 



COS ■ 



ai9). 



4.2. Examples of distributed loads. To check the convergence of the expressions in ()24p we 
will consider a symmetric constant distributed load o"(0) of the form 

^(<^) = ^— y^, O<0<</)o, (36) 



sm 
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-N=4 
-N=10 
-N=100 



N=4 
N=10 
N=100 




(degrees) 



Figure 8. Convergence of the expressions for Ur = ^^^^^ Ur and Uq — — yr 
given in equation (p^ with N for a Hertzian-type distributed load given by ([37|) 
The load was distributed up to = 10°. 



U0 



and a symmetric Hertzian-type load of the form 

a{<l>) 



-3F 



sm^ (/)o 



sm 



< (/) < 



sm 



(37) 



Both loads have been normalized such that their resultant forces are —F for all ranges of the 
angle (j)Q, which is equivalent to the point force given by equation ([7]). The solution on the 
interval < < (f>Q is obtained using (f35l) and for 0o < ^ < vr we apply equations (f29l) directly. 

Firstly, the convergence of the proposed solution (eq. ([35]) ) is compared to the series solution 
for a Hertzian-type load in Figures [5] and [6l These curves indicate that the convergence of the 
radial displacement Ur in the proposed solution is substantially superior to the series solution. 
Figures [7] and [8] show the convergence of the displacements with the truncation limit N under 
both types of loading. Subsequently, Figures [9] and [10] demonstrate that in the limit as (po ^ 
the displacements due to the distributed loads approach those obtained for the point load. 
Moreover, the normalized radial displacement, Ur, is almost indistinguishable from the point 
load for a cpo as large as 10 degrees. Poisson's ratio of z^=0.4 has been used throughout. 

We would also like to investigate how the displacement due to a Hertzian-type load compares 
with that from the Hertzian contact theory. The dimensionless vertical displacement that we 
obtain by the methods outlined in this paper has the form 



U, = Ur COS e-Uesme = (SttGR) 



(38) 



where Uz is the physical vertical displacement. 

The Hertz contact theory [10] is formulated in terms of the radius of the contact area a, the 
displacements directly under the load 5, and the magnitude of the applied load F. We need to 
reformulate these quantities in terms of the contact angle (pQ. The radius of the contact area is 
simply 

a = RsiiKpQ. (39) 
The maximum vertical displacement is related to a in the following manner 

6 = ^ = Rsm^(l)o = 2uM, (^0) 
K 
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Figure 9. Comparison of the displacement due to a constant distributed load 
to the displacement due to a point force of the same magnitude (black). The 
solution was truncated at iV=300. The loads were distributed over 0o=lO°(red), 
30° (green), 50° (blue). 
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Figure 10. Convergence of the displacement due to a Hertzian- type distributed 
load to the displacement due to a point force of the same magnitude (black) . The 
solution was truncated at A^=300. The loads were distributed over (j)Q=W° (red) , 
30° (green), 50° (blue). 



where eq. (I39p was used and the last equality arises from the fact that the Hertzian solution 
presented here is for the contact of two spheres hence we need to half the total displacement. 
Furthermore, Hertz contact theory tells us that the resultant force F is proportional to a^, or 
more accurately 

F=-{- h7=o ^ U?^sin^0o- (41) 



3\l-iyJ R 3Vl 

This allows to rewrite equation ()38p for the dimensionless vertical displacement via Hertz contact 
theory, denoted as (0). Substituting equations ()3U|) and ([^T|) into ([55]) yields 

f/f (0) = (8vrGi^) ,,|r\'^% = (42) 
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350 




Figure 11. Comparison of the dimensionless vertical displacement Uz = 
(SttGR)^ as obtained by the methods in this paper for a Hertzian-type load 
and that obtained from Hertz contact theory defined in ()42p . as a function 
of the contact angle cpQ. 



Equation (I42p gives a way to compare the presented solution for the Hertzian-type load to 
the solution from Hertz contact theory. The numerical results are presented in Figure \TT\ which 
compares the vertical displacements (eq. (j38p with eq. (|42p ) as a function of the contact angle 
(/)Q. Note that along with f/^(0) and U^{0) we also plot 2Uz{4>o), which according to Hertz 
theory should be equal to UziO). The normalized difference between the displacements is shown 
in Figure [T2l As expected, the solutions are close for small contact areas and diverge as this 
area increases. The same can be said about the relationship between the displacements Uz{0) 
and 2Uz{(po)- 

Comparing the maximum displacements ^7^(0) with U^{0) does not tell us anything about 
the shape of the contact area for a sphere loaded by a Hertzian-type load. Hertz contact theory 
states that the contact area between two identical spheres is flat, and thus we can describe 
it using R{cos9 — cos^o)- Therefore, we define a function s{9) to determine how close is our 
calculated displacement to the Hertzian solution as 

s{e) = kUz{9) - {cos 6- cos(^o), (43) 

where A; is a constant determined by enforcing s(0) = s{(j)o), which results in 

= TTTHr^TTTZT^l " ^^^-^o) - (cos0 - cos,^o). (44) 

The function s{9) is plotted in Figure [T3] for several angles (j)Q. These results show that the 
contact area is flat for small contact angles, but gains curvature for larger angles. According to 
Hertz theory, for small contact angles <j)Q, the function s{6) behaves as a constant s{6) ~ </'o/2- 
The angles shown in Figure [T3l are too large to see this behaviour, however, at (j)Q = 5° the values 
are close with s{9) = 0.00334 and (j)l/2 = 0.00381. 
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Figure 12. Normalized difference between the displacements in Figm^e [TT] as a 
function of the contact angle cpo . 



5. Conclusions 

A compact Green's function for a sphere is presented which uses the fundamental idea of 
expressing a slowly convergent series with analytical functions and a quickly convergent series. 
The increased speed of convergence is demonstrated for the point force solution, which is also 
shown to be consistent with the more general distributed loading in the limit as the contact angle 
approaches zero. Since the general Green's function contains elliptical integrals, an easy method 
for dealing with the singularity in the integrand is presented. Comparing the exact displacement 
due to a Hertzian-type distributed load to the displacement given by Hertz contact theory we 
conclude that the Hertz contact theory gives accurate results for contact angles up to about 
10°, with a steadily increasing error. For larger contact angles. Hertz theory overestimates the 
displacements and cannot account for the shape of the contact area. This is to say that the 
stress distribution assumed in Hertz theory results in a curved contact surface for larger contact 
angles. 



Appendix A: Legendre polynomial formulas 
The orthogonality and completeness relations for the Legendre functions are 



71 

i(2n + 1) j Pm{cose)Pnicose) smO dO = 6mn, (45a) 


y;i(2n + l)P„(cos0)P„(cos(/>) = ^^^rf , (45b) 

^ — ^ / Sin 



smc 
n=0 
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Figure 13. The flatness of the area under the load as defined in eq. (j44p as 
a function of the normalized angle /(j)Q. Each curve corresponds to a different 
contact angle ranging from (pQ = 5° to 40° in 5° increments. 



Starting with the definition for Pn{x), 

1 r 

Pn{cos9) = — I (cos + i sin cos a)" da, (46) 
71" Jo 

oo 

and using = \z\ < 1, the well known generating function follows 



n=0 

oo -. 

y^rPJ COS 9)= , |t| < 1. (47) 

^ ' Vl + t2-2icos^ ' ' " ^ ^ 

n=0 

Integrating the identity (j47p with respect to t implies 

oo 



r,4--\ Pn{cos6) . , 1, . n 1 /'t — cos9\ 

\ ^n+i_n\ J_ ^ sinh-i(cot 9) + sinh~^ 

^ n + 1 ^ ^ V sin6l / 

n=0 ^ ^ 



, COS \ 

= log(l + csc-) -sinh-Vtan-) + sinh"M — ), |t| < 1, < 6^ < vr. (48) 

^2' 2 \ sin^ / 

Taking the limit as t — ?• 1 yields (jl2ap . Si{9) of \12h\ follows from a similar result jlH eq. 
5.10.1.4], while S2{9) of ()12cp follows from the recurrence relation 

{n + l)Pn+i{x) - (2n + l)xPn{x) + nP„_i(x) = 0, (49) 

after dividing by n and summing from n = 1 to cxo («S'2 agrees with |1H eq. 5.10.1.6]). The 
recurrence relation can be used to then find Sj{9) for j = 3, 4, . . .. 
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Series of products of Legendre functions given by eqs. 6.11.3.1 and 6.11.3.2 of [l2] 

(50) 

n=l 



for— l<x<y<l. Equation (jSOjl i can be derived by operating on both sides by the 
Legendre differential operator L.^ = ^(1 — ^^)^; a^id using the eigenvalue property LxPn{x) = 
— n{n + l)Pn{x) to arrive at ()8bp (for x < y). At the same time, the constants in the right 
member of (I50p i follow by considering the formula for x = 0, y = 1 in which case the sum on 
the left can be found. Equation (f50|) i gives Sq{9, (j)) + Si{0, (f)) by noting that ^ + = 
The following is a simple consequence of Legendre's addition formula |13l eq. 3.19] 

1 r 

Pn{cos6)Pn{cos(j)) = — / P„ (cos 6* COS — sin sin (/) COS a) da. (51) 

Jo 

Multiply both sides of ([5T]) by and sum, implies, using eq. (jT7|) . the identity [HI eq. 5.10.2.1] 

for |t| < 1, 

oo . 

PJ COS 9) PJ COS 4>) = — -k( '^^^'^~ ),u± = Jl-2tcos(e±(f>)+t^. (52) 

^-^ 7r(u+ + n_) \u++U-/ 

n=0 

Appendix B: Analytical functions and their derivatives 
We require the derivatives with respect to 9 of the functions defined in eq. ()12p . They are 

dSoi9) _ sin f - 1 

d9 ~ sin 9 ' 

dSii9) _ dSo{9) 

d9 d9~ """2 

dS2{9) dS,{9) 9( 9 

— — — = — — — cos 6* - i>i (6^) sm + COS - 2 sm - - 1 
du do 2 \ 2 

Similarly, the analytical function used to find ug (eq. (j29p ) in section U] is 

S„,(9,*)=S„«-.« + S,(9,*) = |-^'°«^='"|~i- (54) 

I — 21ogcos I sm |, 9 < (p. 



cot-, (53) 



The derivative of 5oi(^, (j)) is 



dSm{9,(t)) _ J -cotf, 9> 
89 ~|tan|, 9< 
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